Calculate species vulnerability from traits
Compiled on Sat Oct 8 18:33:44 2022 by ohara
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(oharac) ### remotes::install_github('oharac/oharac')
library(data.table)
library(tidyverse)
library(here)
source(here('common_fxns.R'))
options(dplyr.summarise.inform = FALSE) 1 Summary
Read in taxonomic traits filled in by taxon experts and cleaned in prior scripts, then downfilled and gapfilled. Combine with coded sensitivity, adaptive capacity, and exposure from stressor-trait sheets to calculate vulnerability.
2 Data
_raw_data/xlsx/spp_traits_all.xlsx is the raw workbook
prepared by Nathalie Butt from the various submissions of the taxa-group
experts. This has been processed and cleaned by multiple scripts to
files:
Mazu:spp_vuln/spp_vuln_framework/2_gapfill_traits/spp_match_downfill_spp_traits.csvMazu:spp_vuln/spp_vuln_framework/2_gapfill_traits/spp_match_downfill_levels.csvMazu:spp_vuln/spp_vuln_framework/2_gapfill_traits/spp_match_downfill_traits.csvMazu:spp_vuln/spp_vuln_framework/2_gapfill_traits/spp_match_downfill_taxa.csv
See earlier scripts in the process for details
trait_stressor_rankings/stressors_traits_scored.xlsx is
a workbook with each sheet indicating sensitivity or adaptive capacity;
columns in each sheet indicate stressors, and rows indicate traits.
3 Methods
3.1 Read in cleaned and filled trait values.
Set up a function to consistently clean trait values. Trait values in the species trait file are already cleaned and adjusted in many cases to get around mismatches; they are generally lower case, no punctuation except for greater/less than signs.
This function also cleans up category and trait names for consistency. All lower case, punctuation and spaces replaced with underscores. The species trait file is already cleaned in this manner.
clean_traitnames <- function(df, overwrite_clean_col = FALSE) {
df <- df %>%
mutate(category = str_replace_all(category, '[^A-Za-z0-9]+', '_') %>% tolower(),
category = str_replace_all(category, '^_|_$', ''),
trait = str_replace_all(trait, '[^A-Za-z0-9]+', '_') %>% tolower(),
trait = str_replace_all(trait, '^_|_$', ''))
if(!overwrite_clean_col & ('trait_value' %in% names(df))) {
return(df) ### without overwriting existing trait_value
}
if(overwrite_clean_col & ('trait_value' %in% names(df))) {
x <- readline(prompt = 'Overwriting existing trait_value column? y/n ')
if(str_detect(x, '^n')) stop('dammit!')
}
### overwrite existing, or add new
df <- df %>%
mutate(trait_value = str_replace_all(tolower(trait_value), '[^0-9a-z<>]', ''))
return(df)
}
clean_traitvals <- function(df) {
x <- df$trait_value
### First: remove numeric commas
y <- str_replace_all(x, '(?<=[0-9]),(?=[0-9])', '') %>%
### then: drop all non-alphanumeric and a few key punctuation:
str_replace_all('[^0-9a-zA-Z<>,;\\-\\.\\(\\)/ ]', '') %>%
### lower case; do it after dropping any weird non-ascii characters:
tolower() %>%
str_trim() %>%
str_replace_all('n/a', 'na') %>%
### convert remaining commas and slashes to semicolons:
str_replace_all('[,/]', ';') %>%
### drop spaces after numbers e.g. 3 mm -> 3mm:
str_replace_all('(?<=[0-9]) ', '') %>%
### drop spaces before or after punctuation (non-alphanumeric):
str_replace_all(' (?=[^a-z0-9\\(])|(?<=[^a-z0-9\\)]) ', '') %>%
### manually fix some valid slashes:
str_replace_all('nearly sessile;sedentary', 'nearly sessile/sedentary') %>%
str_replace_all('live birth;egg care', 'live birth/egg care') %>%
str_replace_all('chitin;caco3mix', 'chitin/caco3 mix') %>%
str_replace_all('0.5-49mm', '0.5mm-49mm')
df$trait_value <- y
return(df)
}
assign_rank_scores <- function(x) {
y <- tolower(as.character(x))
z <- case_when(!is.na(as.numeric(x)) ~ as.numeric(x),
str_detect(y, '^na') ~ 0.00,
str_detect(y, '^n') ~ 0.00, ### none, NA, no
str_detect(y, '^lo') ~ 0.33,
str_detect(y, '^med') ~ 0.67,
str_detect(y, '^hi') ~ 1.00,
str_detect(y, '^y') ~ 1.00, ### yes
TRUE ~ NA_real_) ### basically NA
return(z)
}Since the species trait file is already cleaned, DO NOT use the
clean_traitvals function - it will overwrite the
trait_value column. Here we will drop plants and algae as
physiologies are so fundamentally different from animals.
spp_traits <- get_spp_traits() ### from common_fxns.R
# cp <- spp_traits %>% filter(species == 'caulophryne polynema')
# cs <- spp_traits %>% filter(species == 'centrophryne spinulosa')str_trait_f <- here('_raw_data/xlsx',
'stressors_traits_scored.xlsx')
str_trait_shts <- readxl::excel_sheets(str_trait_f)3.2 Calculate sensitivity scores
3.2.1 Determine habitats for habitat loss/degradation stressor
Habitat loss and degradation can be considered as an exposure variable, in the same way as potential exposure above. However, in this case, we consider only one stressor.
Questions to consider that will affect scoring/weighting:
3.2.1.1 Is there an actionable difference between across-stage and within-stage dependence?
- Multiple habitats in the “within-stage” category seems to suggest
that a species can move among habitats therefore being less sensitive to
degradation of one habitat in its range. This is a “parallel habitats”
interpretation.
- However, some species may depend on various habitats in a “series habitat” interpretation, e.g., birds that as adults depend on one habitat type for nesting/breeding, another type for forage, and a third for stopovers in migration. In this case, harm to any would present a bottleneck.
- Multiple habitats in the “across-stage” category seem to indicate a
“series” interpretation - e.g., a fish species whose larvae grow in
mangroves, then adults move to reefs.
- However, this could also indicate stages that could survive in multiple habitats (e.g., parallel).
- Because the trait category is not well defined, we cannot systematically distinguish between series and parallel interpretations for either across- or within-stage dependence.
- A series interpretation would sum the vulnerabilities; a parallel interpretation would take an average. Which is most conservative? Parallel has the advantage that it also avoids overweighting based on the number of habitats scored, but would communicate the less alarming results.
To score this we will simply lump together all unique listed habitats, for both within-stage and for across-stage. Sensitivity to habitat degradation or loss will be based on whether the species has any within-stage dependencies and/or any across-stage dependencies, regardless of which habitats or how many.
3.2.2 Calculate for all
sens_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'sensitivity')
sens_traits_df <- sens_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, sens_score, -category, -trait, -trait_value) %>%
mutate(sens_score_orig = as.character(sens_score),
sens_score = assign_rank_scores(sens_score)) %>%
clean_traitnames() %>%
clean_traitvals() %>%
filter(!is.na(category))
### write out traits that increase sensitivity
x <- sens_traits_df %>%
filter(sens_score > 0 & !is.na(sens_score)) %>%
arrange(stressor)
write_csv(x, 'sens_traits_nonzero.csv')
str_sens_trait_scores <- sens_traits_df %>%
select(category, trait, trait_value, sens_score, stressor) %>%
mutate(sens_score = ifelse(is.na(sens_score), 0, sens_score)) %>%
filter(!is.na(trait))To score for a species/stressor combo, first resolve multiple mutually exclusive trait values (using trait_prob) then sum across all traits.
Fix the habitats - if any dependent habitats, set prob to sum of prob across all habitats and trait value to “habitat list” so it will join. Break habs out into a new column for reference.
3.2.2.1 NOTE: if this dramatically messes up habitat loss vulnerability, just set trait_prob = 1.
spp_traits_hab_fixed <- spp_traits %>%
filter(str_detect(trait, 'dependent_habitat')) %>%
group_by(taxon, species, category, trait) %>%
summarize(dep_habs = paste(trait_value, collapse = ';'),
trait_value = 'habitat list',
trait_prob = sum(trait_prob),
.groups = 'drop') %>%
mutate(trait_prob = trait_prob / max(trait_prob)) %>%
bind_rows(spp_traits %>% filter(!str_detect(trait, 'dependent_habitat')))3.2.3 Process species sensitivity
Break down by taxon and process in a loop to reduce computational pressures…
sens_int_file <- here_anx('int/sensitivity_by_spp_traits.csv')
# if(!file.exists(sens_int_file)) {
### unlink(sens_int_file)
spp_sens_raw <- str_sens_trait_scores %>%
left_join(spp_traits_hab_fixed, by = c('category', 'trait', 'trait_value')) %>%
filter(!is.na(stressor) & !is.na(taxon))
message('Summarizing sensitivity trait scores by spp/stressor/trait...')
summarize_sens_by_taxon <- function(t) {
### spp <- spp_vec[1]
t_df <- spp_sens_raw %>%
filter(taxon == t)
message('Processing ', t_df$species %>% n_distinct(), ' species in taxon ', t, ', sum by trait...')
int_df <- t_df %>%
group_by(species, stressor, taxon, trait) %>%
summarize(sens_score = sum(sens_score * trait_prob, na.rm = TRUE))
message('... and now summarizing sensitivity for ', t_df$species %>% n_distinct(),
' species in taxon ', t, '...')
out_df <- int_df %>%
group_by(species, stressor, taxon) %>%
summarize(sens_score = sum(sens_score, na.rm = TRUE), .groups = 'drop')
return(out_df)
}
taxon_vec <- spp_sens_raw$taxon %>% unique() %>% sort()
spp_sens_list <- parallel::mclapply(taxon_vec, mc.cores = 6, FUN = summarize_sens_by_taxon)
spp_sens <- rbindlist(spp_sens_list)
write_csv(spp_sens, sens_int_file)
# }
spp_sens <- read_csv(sens_int_file)3.2.4 Check matching
Unmatched traits between sensitivity scoring sheets and species trait sheets:
x <- spp_traits %>% select(category, trait, trait_value) %>% distinct()
y <- str_sens_trait_scores %>% select(category, trait, trait_value) %>% distinct()These traits are in the species-trait scoring sheets but not found in the sensitivity trait scores (should be adaptive capacity/exposure traits only):
number_of_sites, sub_population_dependence_on_particular_sites, number_of_sites_incl_terrestrial_wetlands, adult_mobility, planktonic_larval_duration_pld_exposure, thermal_sensitivity_to_ocean_warming_max_temps_tolerated, age_to_1st_reproduction_generation_time, are_there_sub_populations, can_the_sex_ratio_be_altered_by_temperature, fecundity, global_population_size, lifetime_reproductive_opportunities, max_age, parental_investment, post_birth_hatching_parental_dependence, reproductive_strategy, depth_min_max, eoo_range, zone
These traits are in the trait-sensitivity scoring sheet but not found in the species scoring (need to be scored for species):
3.2.5 Sensitivity by species group and stressor
Here sensitivity for each stressor is normalized by max observed for that stressor. Biomass removal is excluded, since it is automatically 1 for every species.
plot_df <- spp_sens %>%
group_by(stressor) %>%
mutate(sens_norm = sens_score / max(sens_score),
sens_norm = ifelse(is.nan(sens_norm), 0, sens_norm))
p <- ggplot(plot_df, aes(x = stressor, y = sens_norm)) +
geom_jitter(size = .25, alpha = .6, height = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5, size = 6)) +
ylim(c(0, 1)) +
facet_wrap(~ taxon)
ggsave(here('figs/spp_sens_scores.png'), height = 6, width = 6, dpi = 300)
knitr::include_graphics(here('figs/spp_sens_scores.png'))
3.2.6 Sensitivity to top three stressors by taxon
Here sensitivity for each stressor is normalized by max observed for that stressor. Biomass removal is excluded, since it is automatically 1 for every species.
top_3_sens <- spp_sens %>%
group_by(stressor) %>%
mutate(sens_norm = sens_score / max(sens_score)) %>%
group_by(taxon, stressor) %>%
summarize(sens_norm = mean(sens_norm) %>% round(3)) %>%
arrange(desc(sens_norm)) %>%
group_by(taxon) %>%
filter(sens_norm >= nth(sens_norm, 3)) %>%
filter(!stressor %in% c('biomass_removal', 'sst_rise'))
DT::datatable(top_3_sens)3.3 Score general adaptive capacity
General adaptive capacity traits are basically related to the overall population’s resilience in the face of a threat. Large extents of occurrence, large population sizes, presence of multiple subpopulations, and reproductive strategies fall into this category.
adcap_gen_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'gen_adcap')
adcap_gen_traits <- adcap_gen_traits_raw %>%
select(category, trait, trait_value, adcap_score) %>%
# filter(trait != 'max age' & trait != 'if one/few, size') %>%
mutate(adcap_score_orig = as.character(adcap_score),
adcap_score = assign_rank_scores(adcap_score)) %>%
clean_traitnames() %>%
clean_traitvals()
spp_adcap_gen_raw <- spp_traits %>%
inner_join(adcap_gen_traits, by = c('category', 'trait', 'trait_value')) %>%
mutate(adcap_gen_score = ifelse(is.na(adcap_score), 0, adcap_score)) %>%
select(-adcap_score, -adcap_score_orig)
spp_adcap_gen <- spp_adcap_gen_raw %>%
mutate(wt_adcap_gen = trait_prob * adcap_gen_score) %>%
group_by(species, taxon) %>%
summarize(adcap_gen_score = sum(wt_adcap_gen, na.rm = TRUE),
.groups = 'drop')3.3.1 General adaptive capacity by taxon
Here scores are not normalized. The pattern would be identical if normalized, simply rescaled.
p <- ggplot(spp_adcap_gen, aes(x = taxon, y = adcap_gen_score)) +
geom_jitter(size = .25, alpha = .6, width = .2, height = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))
ggsave(here('figs/spp_adcap_gen_scores.png'), height = 4, width = 4, dpi = 300)
knitr::include_graphics(here('figs/spp_adcap_gen_scores.png'))
- Median: 7.67
- Mean: 7.8639491
- Standard Deviation: 2.2155742
3.4 Score specific adaptive capacity
Specific adaptive capacity traits are basically related to an organism’s ability to avoid or mitigate exposure, primarily through movement and larval dispersal.
adcap_spec_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'spec_adcap') %>%
janitor::clean_names() %>%
filter(!str_detect(tolower(category), 'spatial')) ### drop exposure traits
### now, clean up the result and assign scores
adcap_spec_traits <- adcap_spec_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, adcap_score, -category, -trait, -trait_value) %>%
mutate(adcap_score_orig = as.character(adcap_score),
adcap_score = assign_rank_scores(adcap_score)) %>%
clean_traitnames() %>%
clean_traitvals() %>%
filter(!is.na(trait))
spp_adcap_spec_raw <- spp_traits %>%
inner_join(adcap_spec_traits, by = c('category', 'trait', 'trait_value')) %>%
mutate(adcap_spec_score = ifelse(is.na(adcap_score), 0, adcap_score)) %>%
filter(!is.na(stressor)) %>%
select(-adcap_score, -adcap_score_orig)
spp_adcap_spec <- spp_adcap_spec_raw %>%
mutate(wt_adcap_spec = adcap_spec_score * trait_prob) %>%
group_by(species, stressor, taxon) %>%
summarize(adcap_spec_score = sum(adcap_spec_score, na.rm = TRUE),
.groups = 'drop')
adcap_spec_sum <- spp_adcap_spec %>%
group_by(stressor) %>%
summarize(median = median(adcap_spec_score, na.rm = TRUE),
mean = mean(adcap_spec_score, na.rm = TRUE),
sd = sd(adcap_spec_score, na.rm = TRUE),
.groups = 'drop')3.4.1 Check matching
Unmatched traits between specific adaptive capacity scoring sheet and species trait sheets:
x <- spp_traits %>% select(category, trait, trait_value) %>% distinct()
y <- adcap_spec_traits %>% select(category, trait, trait_value) %>% distinct()Traits in species-trait sheets, not in specific adaptive capacity scores:
biomass_removal_sensitivity_trait, adult_body_mass_body_size, biomineral, calcium_carbonate_structure_location, calcium_carbonate_structure_stages, communication_requirement_sound, extreme_pressure_wave_sensitive_structures, flight, respiration_structures, number_of_sites, sub_population_dependence_on_particular_sites, number_of_sites_incl_terrestrial_wetlands, dissolved_oxygen, ph, salinity, sensitivity_to_wave_energy_physical_forcing, sst_rise_sensitivity_trait, thermal_sensitivity_to_heat_spikes_heat_waves, thermal_sensitivity_to_ocean_warming_max_temps_tolerated, thermal_tolerance_range, age_to_1st_reproduction_generation_time, are_there_sub_populations, can_the_sex_ratio_be_altered_by_temperature, fecundity, feeding_larva_post_hatching_metamorphosis, global_population_size, lifetime_reproductive_opportunities, max_age, parental_investment, post_birth_hatching_parental_dependence, reproductive_strategy, depth_min_max, eoo_range, zone, across_stage_dependent_habitats_condition, air_sea_interface, dependent_interspecific_interactions, extreme_diet_specialization, habitat_forming, photosynthetic, terrestrial_and_marine_life_stages, within_stage_dependent_habitats_condition, navigation_requirements_light, navigation_requirements_magnetic, navigation_requirements_sound
Traits in specific ad cap scores, not in spp-traits:
can_the_sex_ratio_be_affected_by_temperature
3.4.2 specific adaptive capacity by stressor and species group
Here these are normalized for each stressor.
plot_df <- spp_adcap_spec %>%
group_by(stressor) %>%
mutate(adcap_spec_norm = adcap_spec_score / max(adcap_spec_score),
adcap_spec_norm = ifelse(is.nan(adcap_spec_norm), 0, adcap_spec_norm))
p <- ggplot(plot_df, aes(x = stressor, y = adcap_spec_norm)) +
geom_jitter(size = .25, alpha = .6, width = .2, height = 0) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5, size = 6)) +
ylim(c(0, 1)) +
facet_wrap( ~ taxon)
ggsave(here('figs/spp_adcap_spec_scores.png'), height = 6, width = 6, dpi = 300)
knitr::include_graphics(here('figs/spp_adcap_spec_scores.png'))
| stressor | median | mean | sd |
|---|---|---|---|
| air_temp | 2.01 | 2.6358826 | 1.6824621 |
| biomass_removal | 0.00 | 0.0000000 | 0.0000000 |
| bycatch | 0.00 | 0.0000000 | 0.0000000 |
| entanglement_macroplastic | 1.00 | 1.4119970 | 0.7198893 |
| eutrophication_nutrient_pollution | 1.67 | 2.4591960 | 1.6115212 |
| habitat_loss_degradation | 1.00 | 1.2429333 | 1.0118198 |
| inorganic_pollution | 1.67 | 2.3596339 | 1.4939622 |
| invasive_species | 0.33 | 0.7403439 | 0.8405954 |
| light_pollution | 1.00 | 1.3424955 | 1.1447910 |
| marine_heat_waves | 0.33 | 0.7091713 | 0.7747627 |
| noise_pollution | 0.67 | 1.0949515 | 1.0580643 |
| oa | 1.67 | 2.4591960 | 1.6115212 |
| oceanographic | 0.00 | 0.4235647 | 0.6498276 |
| organic_pollution | 2.00 | 2.7067401 | 1.6911842 |
| plastic_pollution_microplastic | 1.67 | 2.1795405 | 1.3069658 |
| poisons_toxins | 1.67 | 2.4591960 | 1.6115212 |
| salinity | 1.67 | 2.5558299 | 1.7312211 |
| sedimentation | 1.67 | 2.4591960 | 1.6115212 |
| slr | 1.00 | 1.4391294 | 1.2785403 |
| sst_rise | 0.00 | 0.0000000 | 0.0000000 |
| storm_disturbance | 2.01 | 2.7666173 | 1.8483572 |
| uv | 1.34 | 1.8678620 | 1.2423140 |
| wildlife_strike | 0.67 | 1.0949515 | 1.0580643 |
3.4.3 Adaptive capacity to top three stressors by taxon
top_3_adcap <- spp_adcap_spec %>%
group_by(taxon, stressor) %>%
summarize(adcap_score = mean(adcap_spec_score),
.groups = 'drop') %>%
arrange(desc(adcap_score)) %>%
group_by(taxon) %>%
filter(adcap_score >= nth(adcap_score, 3)) %>%
ungroup()
DT::datatable(top_3_adcap)3.5 Assign exposure potential modifier
Exposure potential modifier checks whether the depth and oceanic zones of the stressor match with the depth and oceanic zones of the species. These fall into the “spatial scale” category with the exception of EOO.
exp_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'spec_adcap') %>%
filter(str_detect(tolower(category), 'spatial')) ### include only exposure traits
exp_traits <- exp_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, exp_score, -category, -trait, -trait_value) %>%
mutate(exp_score_orig = as.character(exp_score),
exp_score = assign_rank_scores(exp_score)) %>%
clean_traitnames()
### wildlife strike calculated separately to avoid spurious exposure of
### benthic/intertidal critters
spp_exposure_raw <- spp_traits %>%
inner_join(exp_traits, by = c('category', 'trait', 'trait_value')) %>%
filter(stressor != 'wildlife_strike') %>%
group_by(species, stressor, taxon) %>%
summarize(exposure_mod = as.integer(sum(exp_score, na.rm = TRUE) > 0),
.groups = 'drop') %>%
arrange(stressor, taxon)
air_sea <- spp_traits %>%
filter(trait == 'air_sea_interface' & trait_value != 'no') %>%
.$species %>% unique()
spp_exposure_wildlife_strike <- spp_traits %>%
inner_join(exp_traits, by = c('category', 'trait', 'trait_value')) %>%
filter(stressor == 'wildlife_strike') %>%
group_by(species, taxon) %>%
summarize(exposure_mod = as.integer(sum(exp_score, na.rm = TRUE) > 0),
zone = paste(trait_value, collapse = ';'),
.groups = 'drop') %>%
mutate(btm_only = str_detect(zone, 'benthic|demersal') &
!str_detect(zone, 'oceanic|air') &
!species %in% air_sea,
intertidal_only = str_detect(zone, 'intertidal') & !str_detect(zone, 'neritic|oceanic'),
taxon_excl = taxon %in% c('crustacea_arthropods', 'molluscs'),
clear = btm_only|intertidal_only|taxon_excl) %>%
mutate(override = exposure_mod > 0 & clear,
exposure_mod = ifelse(clear, 0, exposure_mod),
stressor = 'wildlife_strike')
spp_exposure <- spp_exposure_wildlife_strike %>%
select(species, taxon, exposure_mod, stressor) %>%
bind_rows(spp_exposure_raw)
non_exposures <- spp_exposure %>%
group_by(taxon, stressor) %>%
mutate(n_gps = n_distinct(species)) %>%
filter(exposure_mod == 0) %>%
summarize(n_gps_no_exp = n_distinct(species),
n_gps = first(n_gps),
pct_no_exp = round(n_gps_no_exp / n_gps, 3),
.groups = 'drop') %>%
arrange(stressor, desc(n_gps_no_exp))
null_exposures <- spp_traits %>%
select(species, taxon) %>%
anti_join(spp_exposure, by = c('species', 'taxon')) %>%
distinct()
cleared_exposure <- spp_exposure_wildlife_strike %>%
filter(override)
write_csv(cleared_exposure, here_anx('int/cleared_exposure.csv'))3.5.1 These species are not listed as potential exposure to these stressors:
Note: this is exposure potential only, based on overlap between species presence and stressor presence - nothing about sensitivity or actual exposure. Check that these logic out.
3.5.2 These species drop out of the exposure potential calculation
Check the spp traits for these species to identify proper assignment of at least one depth zone or ocean zone.
3.5.3 These species are cleared of otherwise positive wildlife strike
| Var1 | Freq |
|---|---|
| cephalopods | 420 |
| corals | 288 |
| crustacea_arthropods | 6664 |
| echinoderms | 7533 |
| elasmobranchs | 872 |
| fish | 1384 |
| molluscs | 9571 |
| reptiles | 9 |
| sponges | 7436 |
4 Combine scores
We will try a calculation for vulnerability \(V\) of species \(i\) to stressor \(j\) that basically looks like this:
\[\text{sensitivity score } S_{i,j} = \mathbf{s}_j^T \mathbf{t}_i\] based on a vector \(\mathbf{s}_j\) of trait-based sensitivity to stressor \(j\), and vector \(\mathbf{t}_i\) of traits of species \(i\);
\[\text{specific adaptive capacity score } K_{i,j} = \mathbf{k}_j^T \mathbf{t}_i\] based on vector \(\mathbf{k}_j\) of trait-based specific adaptive capacity to stressor \(j\); \[\text{general adaptive capacity score } G_{i} = \mathbf{g}^T \mathbf{t}_i\] based on vector \(\mathbf{g}\) of trait-based general adaptive capacity;
\[\text{exposure potential modifier } E_{i,j} = \begin{cases} 1 \text{ when }\mathbf{e}_j^T \mathbf{t}_i > 0\\ 0 \text{ else} \end{cases}\] based on vector \(\mathbf{e}_j\) of trait-based presence of stressor \(j\) (i.e. depth zones and ocean zones in which stressor occurs).
\[\text{vulnerability } V_{i,j} = \frac{S_{i,j} / {S_j}'}{1 + G_i/ {G}' + K_{i,j}/ {K_j}'} \times E_{i,j}\] Each component (\(S_{i,j}, G_i, K_{i,j}\)) is normalized by a reference value (\(S_{j}', G', K_{j}'\) using mean, median, max, etc) for that component for that stressor across all species. Note: median risks referencing to zero for some stressors with few sensitivities (e.g. light pollution); mean risks having a very low reference for the same. Max risks being driven by an outlier, but here the sensitivity scores are generally capped at some low-ish value since there are a finite number of traits that can confer sensitivity. Therefore, we will use max as the reference point. We may wish to consider max possible, which may differ from max observed, in a future iteration?
For species groups with NA in specific adaptive capacity, force to zero (no matching adaptive traits); for species with NA in exposure potential, force to 1 (assume exposure potential).
These results will be saved by species group for now, for future matching to the species level.
### Check that all stressors are matched to ensure proper combining of scores
exp_strs <- spp_exposure$stressor %>% unique()
sens_strs <- spp_sens$stressor %>% unique()
adcap_strs <- spp_adcap_spec$stressor %>% unique()
if(!all(exp_strs %in% sens_strs) | !all(sens_strs %in% exp_strs)) {
stop('Mismatch between stressors in exposure traits and sensitivity traits!')
}
if(!all(adcap_strs %in% sens_strs) | !all(sens_strs %in% adcap_strs)) {
stop('Mismatch between stressors in ad cap traits and sensitivity traits!')
}
if(!all(exp_strs %in% adcap_strs) | !all(adcap_strs %in% exp_strs)) {
stop('Mismatch between stressors in exposure traits and ad cap traits!')
}4.1 Calc vulnerability
Because our current method imputes traits weighted by probability/distribution, rather than a vulnerability score mean/sd, we can calculate vulnerability scores directly based on traits for each species. This avoids needing to use a Monte Carlo approach as we did for the Ecosphere paper.
Since the vulnerability is based on rescaled sensitivity and ad cap based on the max observed values, we need to identify reference values across the entire set of spp groups.
spp_scores_mean <- spp_sens %>%
left_join(spp_adcap_gen, by = c('taxon', 'species')) %>%
left_join(spp_adcap_spec, by = c('taxon', 'species', 'stressor')) %>%
drop_na()
write_csv(spp_scores_mean, here_anx('int/spp_scores_mean.csv'))
ref_values <- spp_scores_mean %>%
group_by(stressor) %>%
summarize(max_sens = max(sens_score, 1),
### the 1 ensures that stressors with very low scores
### don't get normalized by a low reference
max_adcap_gen = max(adcap_gen_score, 1),
max_adcap_spec = max(adcap_spec_score, 1),
.groups = 'drop')
rescale_vals <- function(x) {
z <- x %>%
left_join(ref_values, by = 'stressor') %>%
mutate(sens_rescale = sens_score / max_sens,
acg_rescale = adcap_gen_score / max_adcap_gen,
acs_rescale = adcap_spec_score / max_adcap_spec)
return(z)
}
calc_vuln <- function(x) {
zz <- x %>%
rescale_vals() %>%
mutate(vuln_raw = sens_rescale / (1 + acg_rescale + acs_rescale))
return(zz)
}spp_vuln_all <- spp_scores_mean %>%
calc_vuln() %>%
select(taxon, species, stressor,
sens_score, adcap_gen_score, adcap_spec_score, vuln_raw) %>%
left_join(spp_exposure, by = c('species', 'taxon', 'stressor')) %>%
mutate(vuln_raw = vuln_raw * exposure_mod) %>%
distinct()
spp_vuln_rescale <- spp_vuln_all %>%
ungroup() %>%
mutate(vuln = vuln_raw / max(vuln_raw)) %>%
arrange(taxon, species, stressor)
spp_vuln_rescale_taxa <- spp_vuln_rescale %>%
select(taxon, species) %>%
distinct() %>%
arrange(taxon, species) %>%
mutate(vuln_tx_id = 1:n())
str_lookup <- spp_vuln_rescale %>%
select(stressor) %>%
distinct() %>%
arrange(stressor) %>%
mutate(vuln_str_id = 1:n())
spp_vuln_rescale_score <- spp_vuln_rescale %>%
left_join(spp_vuln_rescale_taxa, by = c('taxon', 'species')) %>%
left_join(str_lookup, by = 'stressor') %>%
select(vuln_tx_id, vuln_str_id, vuln) %>%
mutate(vuln = round(vuln, 5))Save the vulnerability scores for Github, broken into subtables for file size. This is the info most useful for the typical user.
write_csv(spp_vuln_rescale_taxa, here('_output/spp_vuln_from_traits_tx.csv'))
write_csv(str_lookup, here('_output/spp_vuln_from_traits_str.csv'))
write_csv(spp_vuln_rescale_score, here('_output/spp_vuln_from_traits_score.csv'))
### to reassemble just vuln scores:
# spp_vuln_scores <- fread(here('_output/spp_vuln_from_traits_score.csv')) %>%
# left_join(fread(here('_output/spp_vuln_from_traits_str.csv'))) %>%
# left_join(fread(here('_output/spp_vuln_from_traits_tx.csv'))) %>%
# select(-vuln_tx_id, -vuln_str_id)Write out the full scores (including sensitivity, spec and gen ad capacity, exposure) to server (Mazu). We could save those components too, for convenience for another user, but that user can also just regenerate them with this script…
spp_vuln_resc_all_scores <- spp_vuln_rescale %>%
left_join(spp_vuln_rescale_taxa, by = c('taxon', 'species')) %>%
left_join(str_lookup, by = 'stressor') %>%
select(vuln_tx_id, vuln_str_id,
sens_score, adcap_gen_score, adcap_spec_score,
exposure_mod, vuln)
write_csv(spp_vuln_rescale_taxa,
here_anx('3_vuln_score_traits/spp_vuln_from_traits_tx.csv'))
write_csv(str_lookup,
here_anx('3_vuln_score_traits/spp_vuln_from_traits_str.csv'))
write_csv(spp_vuln_resc_all_scores,
here_anx('3_vuln_score_traits/spp_vuln_from_traits_all_scores.csv'))
### to reassemble all scores incl components:
# spp_vuln_scores_all <- fread(here_anx('3_vuln_score_traits/spp_vuln_from_traits_all_scores.csv')) %>%
# left_join(fread(here_anx('3_vuln_score_traits/spp_vuln_from_traits_str.csv'))) %>%
# left_join(fread(here_anx('3_vuln_score_traits/spp_vuln_from_traits_tx.csv'))) %>%
# select(-vuln_tx_id, -vuln_str_id)
# cp <- spp_vuln_scores_all %>% filter(species == 'caulophryne polynema')
# cs <- spp_vuln_scores_all %>% filter(species == 'centrophryne spinulosa')4.2 Time for some plots!
4.2.1 Vulnerability per stressor by spp gp
strs_to_keep <- c('biomass_removal',
'bycatch',
'eutrophication_nutrient_pollution',
'habitat_loss_degradation',
'light_pollution',
'plastic_pollution_microplastic',
'wildlife_strike',
'marine_heat_waves',
'oa',
'slr',
'sst_rise',
'uv')
plot_df <- spp_vuln_rescale %>%
filter(stressor %in% strs_to_keep) %>%
distinct() %>%
mutate(across(is.numeric, ~round(., 4)))
taxa <- plot_df$taxon %>% unique()
for(t in taxa) { # t <- taxa[6]
t_vuln <- plot_df %>%
filter(taxon == t)
mean_str_vuln <- t_vuln %>%
group_by(stressor) %>%
summarize(vuln = mean(vuln), .groups = 'drop')
mean_tot_vuln <- t_vuln %>%
summarize(vuln = mean(vuln))
vuln_plot <- ggplot(t_vuln,
aes(x = stressor, y = vuln)) +
theme_ohara(base_size = 12) +
geom_hline(data = mean_tot_vuln, aes(yintercept = vuln), color = 'red') +
geom_jitter(size = 1, alpha = .6, width = .2, height = .02) +
geom_point(data = mean_str_vuln,
shape = 21, size = 3,
alpha = 1, color = 'yellow', fill = 'red') +
ylim(0, 1) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
strip.background = element_rect(fill = 'grey90')) +
labs(title = paste0('Vulnerability: ', t))
plotfile <- sprintf('figs/vuln_plot_from_traits_%s.png', t)
ggsave(plot = vuln_plot, filename = plotfile,
width = 8, height = 8, dpi = 300)
# cat(sprintf('\n', plotfile))
}
# knitr::include_graphics(here('figs/vuln_plot.png'))4.2.2 Vulnerability by stressor across all taxa
str_mean_vuln <- plot_df %>%
group_by(stressor) %>%
summarize(vuln = mean(vuln, na.rm = TRUE), .groups = 'drop')
all_mean_vuln <- plot_df %>%
summarize(vuln = mean(vuln, na.rm = TRUE))
x <- ggplot(plot_df, aes(x = stressor, y = vuln, color = taxon)) +
theme_ohara(base_size = 12) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
geom_hline(data = all_mean_vuln, aes(yintercept = vuln), color = 'red') +
geom_jitter(size = 1, alpha = .6, width = .2, height = .02) +
geom_point(data = str_mean_vuln,
shape = 21, size = 3,
alpha = 1, color = 'yellow', fill = 'black') +
ylim(0, 1)
x
4.2.3 Vulnerability per spp gp by stressor
stressors <- plot_df$stressor %>% unique()
for(s in stressors) { # s <- stressors[6]
s_vuln <- plot_df %>%
filter(stressor == s)
mean_str_vuln <- s_vuln %>%
group_by(taxon) %>%
summarize(vuln = mean(vuln), .groups = 'drop')
mean_tot_vuln <- s_vuln %>%
summarize(vuln = mean(vuln))
vuln_plot <- ggplot(s_vuln,
aes(x = taxon, y = vuln)) +
theme_ohara(base_size = 12) +
geom_hline(data = mean_tot_vuln, aes(yintercept = vuln), color = 'red') +
geom_jitter(size = 1, alpha = .6, width = .2, height = .02) +
geom_point(data = mean_str_vuln,
shape = 21, size = 3,
alpha = 1, color = 'yellow', fill = 'red') +
ylim(0, 1) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
strip.background = element_rect(fill = 'grey90')) +
labs(title = paste0('Vulnerability: ', s))
plotfile <- sprintf('figs/vuln_plot_by_stressor_%s.png', s)
ggsave(plot = vuln_plot, filename = plotfile,
width = 8, height = 8, dpi = 300)
# cat(sprintf('\n', plotfile))
}4.2.4 Vulnerability by taxon across all stressors
tx_mean_vuln <- plot_df %>%
group_by(taxon) %>%
summarize(vuln = mean(vuln, na.rm = TRUE), .groups = 'drop')
all_mean_vuln <- plot_df %>%
summarize(vuln = mean(vuln, na.rm = TRUE))
x <- ggplot(plot_df, aes(x = taxon, y = vuln, color = stressor)) +
theme_ohara(base_size = 12) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
geom_hline(data = all_mean_vuln, aes(yintercept = vuln), color = 'red') +
geom_jitter(size = 1, alpha = .6, width = .2, height = .02) +
geom_point(data = tx_mean_vuln,
shape = 21, size = 3,
alpha = 1, color = 'yellow', fill = 'black') +
ylim(0, 1)
x